{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Regression diagnostics"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This example file shows how to use a few of the ``statsmodels`` regression diagnostic tests in a real-life context. You can learn about more tests and find out more information about the tests here on the [Regression Diagnostics page.](https://www.statsmodels.org/stable/diagnostic.html)\n",
    "\n",
    "Note that most of the tests described here only return a tuple of numbers, without any annotation. A full description of outputs is always included in the docstring and in the online ``statsmodels`` documentation. For presentation purposes, we use the ``zip(name,test)`` construct to pretty-print short descriptions in the examples below."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Estimate a regression model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                            OLS Regression Results                            \n",
      "==============================================================================\n",
      "Dep. Variable:                Lottery   R-squared:                       0.348\n",
      "Model:                            OLS   Adj. R-squared:                  0.333\n",
      "Method:                 Least Squares   F-statistic:                     22.20\n",
      "Date:                Tue, 24 Dec 2019   Prob (F-statistic):           1.90e-08\n",
      "Time:                        15:02:23   Log-Likelihood:                -379.82\n",
      "No. Observations:                  86   AIC:                             765.6\n",
      "Df Residuals:                      83   BIC:                             773.0\n",
      "Df Model:                           2                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "===================================================================================\n",
      "                      coef    std err          t      P>|t|      [0.025      0.975]\n",
      "-----------------------------------------------------------------------------------\n",
      "Intercept         246.4341     35.233      6.995      0.000     176.358     316.510\n",
      "Literacy           -0.4889      0.128     -3.832      0.000      -0.743      -0.235\n",
      "np.log(Pop1831)   -31.3114      5.977     -5.239      0.000     -43.199     -19.424\n",
      "==============================================================================\n",
      "Omnibus:                        3.713   Durbin-Watson:                   2.019\n",
      "Prob(Omnibus):                  0.156   Jarque-Bera (JB):                3.394\n",
      "Skew:                          -0.487   Prob(JB):                        0.183\n",
      "Kurtosis:                       3.003   Cond. No.                         702.\n",
      "==============================================================================\n",
      "\n",
      "Warnings:\n",
      "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n"
     ]
    }
   ],
   "source": [
    "from statsmodels.compat import lzip\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import statsmodels.formula.api as smf\n",
    "import statsmodels.stats.api as sms\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Load data\n",
    "url = 'https://raw.githubusercontent.com/vincentarelbundock/Rdatasets/master/csv/HistData/Guerry.csv'\n",
    "dat = pd.read_csv(url)\n",
    "\n",
    "# Fit regression model (using the natural log of one of the regressors)\n",
    "results = smf.ols('Lottery ~ Literacy + np.log(Pop1831)', data=dat).fit()\n",
    "\n",
    "# Inspect the results\n",
    "print(results.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Normality of the residuals"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Jarque-Bera test:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[('Jarque-Bera', 3.393608024843164),\n",
       " ('Chi^2 two-tail prob.', 0.18326831231663396),\n",
       " ('Skew', -0.4865803431122335),\n",
       " ('Kurtosis', 3.003417757881633)]"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "name = ['Jarque-Bera', 'Chi^2 two-tail prob.', 'Skew', 'Kurtosis']\n",
    "test = sms.jarque_bera(results.resid)\n",
    "lzip(name, test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Omni test:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[('Chi^2', 3.713437811597181), ('Two-tail probability', 0.15618424580304824)]"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "name = ['Chi^2', 'Two-tail probability']\n",
    "test = sms.omni_normtest(results.resid)\n",
    "lzip(name, test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Influence tests\n",
    "\n",
    "Once created, an object of class ``OLSInfluence`` holds attributes and methods that allow users to assess the influence of each observation. For example, we can compute and extract the first few rows of DFbetas by:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-0.00301154,  0.00290872,  0.00118179],\n",
       "       [-0.06425662,  0.04043093,  0.06281609],\n",
       "       [ 0.01554894, -0.03556038, -0.00905336],\n",
       "       [ 0.17899858,  0.04098207, -0.18062352],\n",
       "       [ 0.29679073,  0.21249207, -0.3213655 ]])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from statsmodels.stats.outliers_influence import OLSInfluence\n",
    "test_class = OLSInfluence(results)\n",
    "test_class.dfbetas[:5,:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Explore other options by typing ``dir(influence_test)``\n",
    "\n",
    "Useful information on leverage can also be plotted:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeoAAAFnCAYAAABpQwo8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAADrhJREFUeJzt3F+IpXd9x/HPt5sIYZPqhg4rAREWQsUqUTtEt65lEozohUhD1IBEMC2LReJ1pBFBQgtBcqNEXI1FUErjhYFApNs/BKPEllkkUP9AvdggWxZXE7NdL3oRfr3YYzcOM3ueGc8cv5nzet3kmZ3fOfvll8O+93n2mafGGAEAevqD3/cAAMDOhBoAGhNqAGhMqAGgMaEGgMaEGgAaE2oAaGxSqKvqaFU9fZXvX1tVT1TV96rq3sWNBwCrbW6oq+pIkq8lOXyVZfclOTPGeGeSu6rqhgXNBwArbcoZ9UtJPpzk4lXWbCR5bHb8nSTrv9tYAECSXDNvwRjjYpJU1dWWHU5ybnb8fJKjWxdU1ckkJ5Pk8OHDf/qGN7xht7MCwCvWmTNnfjHGWNvt6+aGeqJLSa5L8mKS62df/5Yxxqkkp5JkfX19bG5uLui3BoD+quq5vbxuUXd9n0lyYnZ8S5KzC3pfAFhpuz6jrqrbk7xxjPGFl/3y15I8WVXvSvLGJP++oPkAYKVNPqMeY2zM/vtvWyKdMcZzSe5I8r0k7x5jvLTIIQFgVS3q36gzxvjvXLnzGwBYAE8mA4DGhBoAGhNqAGhMqAGgMaEGgMaEGgAaE2oAaEyoAaAxoQaAxoQaABoTagBoTKgBoDGhBoDGhBoAGhNqAGhMqAGgMaEGgMaEGgAaE2oAaEyoAaAxoQaAxoQaABoTagBoTKgBoDGhBoDGhBoAGhNqAGhMqAGgMaEGgMaEGgAaE2oAaEyoAaAxoQaAxoQaABoTagBoTKgBoDGhBoDGhBoAGhNqAGhMqAGgMaEGgMaEGgAaE2oAaEyoAaAxoQaAxoQaABoTagBoTKgBoDGhBoDGhBoAGhNqAGhMqAGgMaEGgMaEGgAaE2oAaEyoAaCxSaGuqker6pmqemCH7x+pqierarOqvrTYEQFgdc0NdVXdmeTQGON4kmNVdfM2y+5J8o0xxnqSG6pqfcFzAsBKmnJGvZHksdnx6SQntlnzyyRvqqrXJHldkp8tZDoAWHFTQn04ybnZ8fNJjm6z5rtJXp/kk0l+PFv3W6rq5OzS+OaFCxf2OC4ArJYpob6U5LrZ8fU7vOYzST4+xvhskp8k+djWBWOMU2OM9THG+tra2l7nBYCVMiXUZ3LlcvctSc5us+ZIkjdX1aEkb08yFjIdAKy4KaF+PMk9VfVwkg8l+WFVPbhlzd8lOZXkxSQ3JvmHhU4JACvqmnkLxhgXq2ojyR1JHhpjnE/y7JY1/5HkT/ZlQgBYYXNDnSRjjBdy5c5vAGBJPJkMABoTagBoTKgBoDGhBoDGhBoAGhNqAGhMqAGgMaEGgMaEGgAaE2oAaEyoAaAxoQaAxoQaABoTagBoTKgBoDGhBoDGhBoAGhNqAGhMqAGgMaEGgMaEGgAaE2oAaEyoAaAxoQaAxoQaABoTagBoTKgBoDGhBoDGhBoAGhNqAGhMqAGgMaEGgMaEGgAaE2oAaEyoAaAxoQaAxoQaABoTagBoTKgBoDGhBoDGhBoAGhNqAGhMqAGgMaEGgMaEGgAaE2oAaEyoAaAxoQaAxoQaABoTagBoTKgBoDGhBoDGhBoAGhNqAGhMqAGgMaEGgMYmhbqqHq2qZ6rqgTnrHqmq9y9mNABgbqir6s4kh8YYx5Mcq6qbd1j3riSvHWM8seAZAWBlTTmj3kjy2Oz4dJITWxdU1bVJvpzkbFV9YGHTAcCKmxLqw0nOzY6fT3J0mzUfTfKjJA8lubWq7tu6oKpOVtVmVW1euHBhr/MCwEqZEupLSa6bHV+/w2vemuTUGON8kq8nuW3rgjHGqTHG+hhjfW1tba/zAsBKmRLqM7lyufuWJGe3WfPTJMdmx+tJnvudJwMAcs2ENY8nebqqbkryviR3V9WDY4yX3wH+aJKvVtXdSa5NctfiRwWA1TM31GOMi1W1keSOJA/NLm8/u2XN/yT54L5MCAArbMoZdcYYL+TKnd8AwJJ4MhkANCbUANCYUANAY0INAI0JNQA0JtQA0JhQA0BjQg0AjQk1ADQm1ADQmFADQGNCDQCNCTUANCbUANCYUANAY0INAI0JNQA0JtQA0JhQA0BjQg0AjQk1ADQm1ADQmFADQGNCDQCNCTUANCbUANCYUANAY0INAI0JNQA0JtQA0JhQA0BjQg0AjQk1ADQm1ADQmFADQGNCDQCNCTUANCbUANCYUANAY0INAI0JNQA0JtQA0JhQA0BjQg0AjQk1ADQm1ADQmFADQGNCDQCNCTUANCbUANCYUANAY0INAI0JNQA0JtQA0JhQA0BjQg0AjU0KdVU9WlXPVNUDc9YdraofLGY0AGBuqKvqziSHxhjHkxyrqpuvsvxzSa5b1HAAsOqmnFFvJHlsdnw6yYntFlXV7Ul+neT8QiYDACaF+nCSc7Pj55Mc3bqgql6V5NNJ7t/pTarqZFVtVtXmhQsX9jIrAKycKaG+lCuXs6/f4TX3J3lkjPGrnd5kjHFqjLE+xlhfW1vb/aQAsIKmhPpMrlzuviXJ2W3WvDvJJ6rqqSRvqaqvLGQ6AFhx10xY83iSp6vqpiTvS3J3VT04xvj/O8DHGH/+m+OqemqM8VeLHxUAVs/cUI8xLlbVRpI7kjw0xjif5NmrrN9Y2HQAsOKmnFFnjPFCrtz5DQAsiSeTAUBjQg0AjQk1ADQm1ADQmFADQGNCDQCNCTUANCbUANCYUANAY0INAI0JNQA0JtQA0JhQA0BjQg0AjQk1ADQm1ADQmFADQGNCDQCNCTUANCbUANCYUANAY0INAI0JNQA0JtQA0JhQA0BjQg0AjQk1ADQm1ADQmFADQGNCDQCNCTUANCbUANCYUANAY0INAI0JNQA0JtQA0JhQA0BjQg0AjQk1ADQm1ADQmFADQGNCDQCNCTUANCbUANCYUANAY0INAI0JNQA0JtQA0JhQA0BjQg0AjQk1ADQm1ADQmFADQGNCDQCNCTUANCbUANCYUANAY5NCXVWPVtUzVfXADt9/dVV9u6pOV9W3qupVix0TAFbT3FBX1Z1JDo0xjic5VlU3b7PsI0keHmO8J8n5JO9d7JgAsJqumbBmI8ljs+PTSU4k+a+XLxhjPPKyL9eS/HwRwwHAqpty6ftwknOz4+eTHN1pYVUdT3JkjPH9bb53sqo2q2rzwoULexoWAFbNlFBfSnLd7Pj6nV5TVTcm+XySe7f7/hjj1BhjfYyxvra2tpdZAWDlTAn1mVy+3J0ktyQ5u3XB7Oaxbyb51BjjuYVNBwArbkqoH09yT1U9nORDSX5YVQ9uWfOXSd6W5G+q6qmq+vCC5wSAlTT3ZrIxxsWq2khyR5KHxhjnkzy7Zc0Xk3xxXyYEgBU25a7vjDFeyJU7vwGAJfFkMgBoTKgBoDGhBoDGhBoAGhNqAGhMqAGgMaEGgMaEGgAaE2oAaEyoAaAxoQaAxoQaABoTagBoTKgBoDGhBoDGhBoAGhNqAGhMqAGgMaEGgMaEGgAaE2oAaEyoAaAxoQaAxoQaABoTagBoTKgBoDGhBoDGhBoAGhNqAGhMqAGgMaEGgMaEGgAaE2oAaEyoAaAxoQaAxoQaABoTagBoTKgBoDGhBoDGhBoAGhNqAGhMqAGgMaEGgMaEGgAaE2oAaEyoAaAxoQaAxoQaABoTagBoTKgBoDGhBoDGhBoAGhNqAGhMqAGgMaEGgMaEGgAamxTqqnq0qp6pqgd+lzUAwO7MDXVV3Znk0BjjeJJjVXXzXtYAALs35Yx6I8ljs+PTSU7scQ0AsEvXTFhzOMm52fHzSd62lzVVdTLJydmX/1tV/7m7UdmlP0ryi9/3ECvAPu8/e7z/7PFy/PFeXjQl1JeSXDc7vj7bn4XPXTPGOJXkVJJU1eYYY33X0zKZPV4O+7z/7PH+s8fLUVWbe3ndlEvfZ3LlUvYtSc7ucQ0AsEtTzqgfT/J0Vd2U5H1J7q6qB8cYD1xlzTsWPyoArJ65Z9RjjIu5fLPY95PcNsZ4dkukt1vz4py3PbWnadkNe7wc9nn/2eP9Z4+XY0/7XGOMRQ8CACyIJ5MBQGP7GmpPNNt/8/avql5dVd+uqtNV9a2qetWyZzwIpn5Oq+poVf1gWXMdJLvY40eq6v3LmusgmfDnxZGqerKqNqvqS8ue76CY/Tnw9FW+f21VPVFV36uqe+e9376F2hPN9t/E/ftIkofHGO9Jcj7Je5c540Gwy8/p53LlRxWZaOoeV9W7krx2jPHEUgc8ACbu8T1JvjH7Ua0bqsqPbO1SVR1J8rVcfr7ITu5LcmaM8c4kd1XVDVd7z/08o96IJ5rtt43M2b8xxiNjjH+efbmW5OfLGe1A2ciEz2lV3Z7k17n8FyJ2ZyNz9riqrk3y5SRnq+oDyxvtwNjI/M/xL5O8qapek+R1SX62nNEOlJeSfDjJxaus2ciV/xffSXLVvxDtZ6i3Pq3s6B7XsLPJ+1dVx5McGWN8fxmDHTBz93n2TwqfTnL/Euc6SKZ8lj+a5EdJHkpya1Xdt6TZDoope/zdJK9P8skkP56tYxfGGBcn/OTTrtq3n6FeyBPNuKpJ+1dVNyb5fJK5/xbCtqbs8/1JHhlj/GppUx0sU/b4rUlOjTHOJ/l6ktuWNNtBMWWPP5Pk42OMzyb5SZKPLWm2VbOr9u1nGD3RbP/N3b/Zmd43k3xqjPHc8kY7UKZ8Tt+d5BNV9VSSt1TVV5Yz2oExZY9/muTY7Hg9ic/z7kzZ4yNJ3lxVh5K8PYmf390fu2rfvv0cdVX9YZKnk/xrZk80S/LBlz8sZZs175hwyYCZiXv810n+Nsmzs1/64hjjH5c96yvZlH3esv6pMcbG8iZ85Zv4Wb4hyVdz+TLhtUnuGmOc2+bt2MbEPb41yd/n8uXvZ5L8xRjj0u9h3Fe83/w5MLt35Y1jjC+87HuvT/Jkkn9J8me53L6Xdnyv/XzgyezutzuSfGd2uWpPa9iZ/VsO+7z/7PH+s8d9zB65fSLJP807QfVkMgBozM1bANCYUANAY0INAI0JNQA0JtQA0JhQA0Bj/weQZe3EQ8c0UwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 576x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAERCAYAAACQIWsgAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3X2clXWd//HXBxxwAHFQCX4zgIo/FrHIaEcUIRlNI9w1iVqoBH52I9hva1VaCsJKk1UMQ93W/HmPtWii4ISL9yHiDSTgSEpBuirGIBuoAyFj0Pj5/XFdA4dzN9cZzjnXmZn38/HgwTnXXDefc52Z63Nd31tzd0RERBJ1ijsAEREpPUoOIiKSQslBRERSKDmIiEgKJQcREUmh5CAiIimUHIrEzK4ws7lxx9HWmNkCM1ua8P5NM6sp4vFXpDuemS0xs9PzeJyitik3sxvNbGKWn79pZscdwv6vMLMrWrt9qQt/Ly+MO45COizuAEQiONfM+rv7n+IOpJm7j487hkPh7pfEHYOUNj05SFvwPvD1uIMQ6UiUHEqAmX3fzN4ys81mdl64bJ6ZfSthnVfNrCrT+uHyN83sLDN7yszuTlh+pZnVh9tMTlg+3czeNrOnzWypmc0Jl3/VzF4Pf3ZRC7H/s5ldl/D+N2Y2IuG4W8P9/N9DOEW/BL5uZp0TjtPJzH4afq71ZnZKws/czD5pZmvN7KqEZb8wszfM7Doz2xGeq25m9mAY4ytmNixKQMnFTWb2WTPbFv5rMLO6hJ+l+367mdn94fm5PuIxM32/ab+vbOc/uVgkPJ83h+suBMrC5Rea2YJ0nzvT71WW+KeZ2Z/M7M/Nv2utWd/MrgrP83Ize8nMasJ/K9J9PjO7OIyx3sy+l7TON83sTjN7NWH5cDOrC497q5lZuPwb4fl8HqhqKf42z931rwj/gCuAuWmWjwV+A3QHBgNvE/xhngosC9f5O+C5bOuHP3sTeBEYBRwRLhsALA/XrwS2hcuPILgjPxK4BpgTLv8o8DJwFNAHqAf6ZPlcfYGXw9c9gP8GLNx+L9ALOBpY3MrztgC4GHgG+Fz4GWuAb4Tn4XDgrHB513AbD392MtAtYdkoYC1wOXAd8EPg88DPw5i/BCxKOv4KoCZNXJmWdw5/NqmF7/e7wIMEN2jfCf4UWzwX6b7ftN9XS+c/PK8XJryfAPw2PJ9fCM/XccCFwILkz53p9yrp9/2KpGW7gI+Fx7i/+TNk+bwp6wOnh79jFcAngaYwnhpgRfLnC7d9LoyxG7Ad6JGwzlsET6VHhcu6AK+GvzuHA0+EvyNVwDvAseF5+Uvi+WuP/1TnEL+zgVMIfuEh+AWudPffmtlAMzuc4ALzQLb1gc3h+5+4+7PNO3f3t8zsUoIL0JkEFw6Av4X/ygj+IP4aLj8LGAj8PnxfTnBR+590wbv7NjN7z8wGAMOA/3J3N7OdwCbgeuBR4P/kdFZS/T9gWsL7scBt7v4BsDw83lCCiz/AbHdfn7SPVcDu8P9PAZ3c/UEz20OQLMYCfz7EOOcAG939P8P3mb6v04GF7v6hmd0RHj+Kg75fMn9fz5Hb+T8deCA8n4vNrCHDegZZf6+yeRb4N6AW+Ka7/yXX9c3sVOARd28AXjSz5O84Oc4PzGwKMJngOz8KOIbg9wDgYXe/I2G7wQQX/8fC912AkwiS5Wp33wzBE3KEz9umqVgpfgb8m7v3dfe+BHdk9eHPHgVGA58FFkdYH2D1QTs3+xTBHerrBHdSzT4kuJCuBaqBGxP2/4uE/fdL3mcaiwkurPuTmLs3EVwUHwg/Q52ZdWlhP9k8APw9wV1wM096vf+9u6fEHMYEwd0mAGY2m+Au/mng+4cQH2Z2LvAZILGyN9P3ZQnxfpjDYZI/V9rvqxXnPzGebDE1F21m+r3K5nPAvxNcgF8xs96tWN/DWJtluoY1x3kCsBJ4lyCRJTdqSHc+X0s4n5XADUQ/P+2GkkP8ngQmmFlPM6vkwCMzBH/YXwDK3f2tCOuncyrwAnAvcG7S8veB4919tLs3PxksB8aaWV8zOwJYT3DnlM1iggT2CYI7Vszs7wiKU34DfI+g+OnoTDtoibv/FbiboOgK4BGCeoiuZjaa4By80opdjyAosniSoPigVcInp58BE8JYm2X6vl4AJppZJ+CrrT0uGb6vVpz/F4Dx4fk8n+AOG4Kinf7hZzwXOCFcnun3Ki0z60bw/bxIUJy3G/jfrVj/eYLWaxVmVk3wtNgcZz8LDCVIiBA8zb4J3EmQZPq1EOpGoJuZfSr8bn5JUIT5InCamfULv+tPt/SZ2zolh+Kabma7E/6Nc/eHCR6bXyG4sH7b3XeE6z8PjAGWNe+ghfXTeYCg3HYrwUV+d3jheJHgTnybmf23mf3KzHq6+yvAVQRFL78HbnL3l7J9KHffQvCovs7dPwyX/ZGgnuAN4I/Af7j722bW38x+F/F8JbuFA3dvdxKUtb9OcHf5T0kX5ah+TnDx+T3QAAy0hIrvHHyV4C7zOTtQMd0ty/f17wRFem9z4AKXs0zfV6bzn2VXvwL+QFAG/zVgW7j8UaBrWNl7LkFRD2T+vcoU5x6Cc/0ywd37MwTJJaf13f0FgpuEjQR1Zc3FSnXhus8RJMPacPmT4f//Q1Cn9AZBHV6m4+4FJhIk+nrgA+DmsDjpR8Aa4NdA1r+J9sDCShjpYMzsX4Ce7j7HzMoI7v7vcPdfxxyaSGRh0rrC3VfEHEq7oyeHjutp4PNmtpXg7vt9giIIERE9OYiISCo9OYiISAolh5hY0Kt1QdxxiIik02aLlY455hg/7rjj4g6jVfbt28eGDRuoqKigrX4GEWmb1q1bt8PdW+pj0nZ7SB933HGsXbu25RVL0KRJkzjppJPo2rUrCxYsiDscEelAzGxzy2upWKnonnjiCbZs2cLUqVPjDkVEJCMlhyJqbGxk+vTp3HbbbYQDPYqIlCQlhyK68sormTx5MoMGDYo7FBGRrNpsnUNb9NBDD7Fjxw7mz5/P3r17aWxsZNeuXSxZsiTu0EREDqLkUEQbNmzY/3rFihUsWLBAFdIiUpJUrCQiIin05BCTmpoaampq4g5DRCQtPTmIiEgKJQcREUlRkORgZneY2SozuzzLOn3M7Jk0yz9mZk8UIi4REYkm73UOZjYe6OzuI8zsTjMb5O6vJq3Ti2A2p+5Jyw2YTzBDVrtUW1fPvMc2sbWhkcqKcmaMGcy4YVVxhyUicpBCPDnUAIvC148Do9Ks00QwFd+upOVfBZ7KtGMzm2pma81s7fbt2/MQanHV1tUza8nL1Dc04kB9QyOzlrxMbV193KGJiBykEMmhO8HcqwDvAn2SV3D3Xe6+M3GZmR0NTAKuy7Rjd7/V3avdvbp37xYHFSw58x7bROO+poOWNe5rYt5jm2KKSEQkvUIkh91Aefi6Rw7HmAvMcvd9BYipJGxtaMxpuYhIXAqRHNZxoCjpZODNiNuNBq4NJwz/hJnNyX9o8aqsKM9puYhIXAqRHGqByWY2H5gAbIhyoXf3v3P3GnevAV5y94wtndqqGWMGU17W+aBl5WWdmTFmcEwRiYikl/fWSu6+y8xqgHOAn7j7NmB9hnVrclne1jW3SlJrJREpdQUZPsPd3+NAi6WSde2113L99dfvf797926++93v8pGPfIQrrriCfv36cf/993P88cfn7ZjjhlUpGYhIyWuzc0hXV1d7vqcJHTVqFDfddBPnnnsua9asYfPmzcyZM4dly5bl9TgiInExs3XuXt3Sehp4L7R8+XJ69+7N0qVLmTJlCpWVlVRWVrJjxw7ef/99unfv3vJORETaCY2tFLrhhhu49NJL2bJlCx//+Mf3L6+srGTz5kjzcYuItBtKDsCWLVt4/fXXGT16NE1NTfTs2XP/z7p3705DQ0OM0YmIFJ+SA3DPPfcwYcIEAHr16nVQMmhsbKRTJ50mEelYdNUDFi9ezLhx4wCorq5m1apVALg7L774IlVVal0kIh1Lh6+Qfuedd3jjjTf21zOMHTuW6dOnc9ZZZ/HKK69w9NFH079//5ijFBEprg6fHFauXMnw4cP3v+/ZsycPPPAAM2bMoGvXrixcuDDG6ERE4qF+DiIiHUjUfg6qcxARkRRKDiIikqLD1zkk0hSeIiIBJYdQ8xSezTO1NU/hCShBiEiHo2KlkKbwFBE5QMkhpCk8RUQOUHIIaQpPEZEDlBxCmsJTROQAVUiHNIWniMgBSg4JNIWniEhAxUoiIpJCyUFERFIUJDmY2R1mtsrMLs+yTh8zeybh/QAzW2Fmy83sVjOzQsQmIiIty3tyMLPxQGd3HwEMNLNBadbpBdwNdE9YPA34prufBfQHhuY7NhERiaYQTw41wKLw9ePAqDTrNAETgV3NC9x9trv/IXx7NLAjeSMzm2pma81s7fbt2/MatIiIHFCI5NAdqA9fvwv0SV7B3Xe5+850G5vZRGCDu29Ns92t7l7t7tW9e/fOZ8wiIpKgEE1ZdwPN3Yp7kEMCMrOBwL8CZxcgLhERiagQTw7rOFCUdDLwZpSNwnqIe4GvZXqqEBGR4ihEcqgFJpvZfGACsMHM5kTYbiYwAPhZ2GppdAFiExGRCAoyh3T4FHAOsNLdt+X9AGgOaRGR1og6h3RBhs9w9/c40GJJRETaGPWQFhGRFEoOIiKSQslBRERSdOjkMGTIEPr27Uvfvn3p16/f/uUNDQ0MGDCA119/PcboRETi02Hnc9izZw9mxrZtqY2pLr/8cr785S8zcODAGCITEYlfh00O69evZ+jQ1LH96urqWLJkCRs3bowhKhGR0tBhi5Xq6up49tlnqays5IQTTmDp0qUAXHLJJfTs2ZNJkyZx8803U4h+ICIipa7DPjn06NGDa6+9lkmTJrF69WrOP/98jjzySJ5//nnmzZtH//79mTlzJmbGxRdfHHe4IiJF1WGTw5QpU/a/Pu200zj22GNZvXo15513HpdddhkA77//PosWLVJyEJEOp8MWK/3iF79g7969+9/X1wejjA8YMGD/ssMPP5yKioqixyYiErcOmxxWrlzJT3/6U/bs2cNNN91Ely5dGDduHI8++ijbt2+nsbGR22+/ndGjNf6fiHQ8HTY5XH311Tz66KP06dOHRYsWsXTpUgYPHszs2bMZOXIk/fr1o6qqiq9//etxhyoiUnQFGZW1GDQqq4hI7qKOytphnxxERCQzJQcREUmh5CAiIik6bD+HZLV19cx7bBNbGxqprChnxpjBjBtWFXdYIiKxUHIgSAyzlrxM474mAOobGpm15GUAJQgR6ZBUrATMe2zT/sTQrHFfE/Me2xRTRCIi8VJyALY2NOa0XESkvVNyACorynNaLiLS3hUkOZjZHWa2yswuz7JOHzN7JuF9mZk9ZGbPmdnXChFXJjPGDKa8rPNBy8rLOjNjzOBihiEiUjLynhzMbDzQ2d1HAAPNbFCadXoBdwPdExZ/G1jn7iOBL5rZEfmOLZNxw6q4ZvxQqirKMaCqopxrxg9VZbSIdFiFaK1UAywKXz8OjAJeTVqnCZgI/Dppu5nh65VANfBU4kZmNhWYCgePnpoP44ZVKRmIiIQKUazUHagPX78L9Elewd13ufvOVmx3q7tXu3t179698xhyqrvuuosLL7wQgCFDhtC3b1/69u1Lv379CnpcEZFSUIjksBtorsntkcMxWrtd3m3fvp0ZM2YAsGfPHsyMbdu2sW3bNrZs2RJXWCIiRVOIC/A6gqIkgJOBNwu8Xd5ddtllnHPOOQCsX7+eoUOHxhWKiEgsCpEcaoHJZjYfmABsMLM5Eba7G7jSzG4ETgJ+W4DYWvTEE0+wZcsWpk6dCkBdXR3PPvsslZWVnHDCCSxdujSOsEREiirvycHddxFULq8GznT39e6etkmru9ckvN4MnAM8B5zt7k3ptimkxsZGpk+fzm233YaZAdCjRw+uvfZatm7dysKFC7nooov461//WuzQRESKqiBjK7n7exxosZTLdltbs12+XHnllUyePJlBgwbtn1N6ypQp+39+2mmnceyxx7J+/XqGDx8eV5giIgWngfcSPPTQQ+zYsYP58+ezd+9eGhsb+cUvfsEHH3xAly5dAKivr6dTJ3UsF5H2TckhwYYNG/a/XrFiBQsWLOCwww7jpz/9KZdccgl33XUXXbp0UQW1iLR7ugVuwdVXX82jjz5Knz59WLRoEUuXLqVr165xhyUiUlDm7nHH0CrV1dW+du3auMMQEWlTzGydu1e3tJ6eHEREJIWSg4iIpFByEBGRFGqtlEFtXT3zHtvE1oZGKivKmTFmsEZtFZEOQ8khjdq6emYteXn/vNL1DY3MWvIygBKEiHQIKlZKY95jm/YnhmaN+5qY99immCISESkuJYc0tjY05rRcRKS9UXJIo7KiPKflIiLtTeTkYGYfM7MxZjbEzHoUMqi4zRgzmPKyzgctKy/rzIwxg2OKSESkuCJVSJvZz4BK4HjgB8C1wOcKGFesmiud1VpJRDqqqK2Vhrp7jZktd/dlZvbdgkZVAsYNq1IyEJEOK2qx0nYz+yHQy8z+D7CtgDGJiEjMoiaHKcBOYBVwJPDVgkUkIiKxi1qsdIm7z21+Y2Ynmdkx7r6yQHGJiEiMoj45DDWz1Wb2pfD9D4B/LVBMIiISs6jJYSAwCvh2+P4jQFlBImqj1qxZw/nnn8+YMWN45JFH4g5HROSQRC1Weg+4CTjczM4H/g7IOJaEmd0BnAQsc/c5UdYxs17AQoLEs87dp0X/GPHauXMnF1xwAbfccgtmxoQJE9i0aRO9evWKOzQRkVaJ+uQwHrgZ+CzQExgL/CrdimY2Hujs7iOAgWY2KOI6k4GF4QxFR5hZizMVlYq3336bOXPmcOaZZ1JTU0P//v3ZsmVL3GGJiLRapCcHd//AzOqBcuBpoMrdb8+weg2wKHz9OEFx1KsR1nkH+JiZVQD9gT8l79jMpgJTAQYMGBAl9KI48cQTOfHEE2lqaqK2tpZ9+/Zx0kknxR2WiEirRXpyCIuA7gN+DdwLzMuyenegPnz9LtAn4jrPAscC/wL8IVx+EHe/1d2r3b26d+/eUUIvqhtvvJFJkyYxbdo0Onfu3PIGIiIJGhoaGDBgAK+//joAr732GhMnTuTTn/40v/zlL4saS9Ripf9NUKT0GjAa+DDLursJnjAAemQ4Rrp1fgRc7O4/BjbSBvtSTJ8+nVdffZWrr76a1157Le5wRKSNufzyy/nyl7/MwIEDaWpqYvz48UycOJF58+bx/e9/n40bNxYtlqjJYQ/waaAz8E9AtprWdQTFRAAnA29GXKcXQZPZzsCpgEeMLXavvfYav/vd7wDo168fp5xyCps2ae4HEYmurq6OJUuWMHv2bAB27NjBt771LcaPH88nP/lJhg8fzquvJpfQF07U5PBFgnqDy4AhwP/Nsm4tMNnM5gMTgA1mltxiKXmdZcA1wK0EPbGPIii+ahO2bt3KxIkTaWho4O2332bNmjUMGzYs7rBEpA255JJL6NmzJ5MmTeLmm2/mIx/5CFOnTsXdefrpp1m/fj1nnHFG0eKJWiH9PkGREsAPW1h3l5nVAOcAP3H3bcD6FtbZCbwAfDSn6EvEGWecwZQpUxgyZAjdu3fnhhtuoLKyMu6wRKSNeP7553n++eeZN28e/fv3Z+bMmZgZF198MYsWLeKiiy5i2rRpHHnkkUWLydxbLr0xs0fcfWwR4omsurra165dG3cYIiKH7Prrr2flypU8+OCDANx9990sWrSIZcuWAUFF9ahRo/iP//gPampqDulYZrYu7DKQVdRipZfDzm8iIpJnPXv2PKh5/uGHH84HH3zAqlWrAKioqODss89mw4YNRYspanI4BfiVmb1gZk+Z2fJCBiUi0pGMGjWKRx99lO3bt9PY2Mjtt9/OhAkTGD9+PFu2bGHXrl08/vjjDB8+vGgxRa1zOLPQgbRltXX1mjVORFpt8ODBzJ49m5EjR/LOO+9w3nnn8Y1vfIOuXbsyYsQIOnXqxIwZMzjllFOKFlOkOgcI5pAGqoC3gD+5++5CBtaSUqlzqK2rZ9aSl2nc17R/WXlZZ64ZP1QJQkRKTl7rHMI5pK8kaG46ELjn0MJrP+Y9tumgxADQuK+JeY+pn4OItF2R53Nw9y8ADe6+jGA2OAG2NjTmtFxEpC2IOmS35pDOoLKinPo0iaCyojzN2iIi2ZVKHabmkD5EM8YMprzs4EH2yss6M2PM4JgiEpG2qrkOs76hEQfqGxqZteRlauvqW9w236I+OZwL3OruKitJ0pzRSyHTi0jblq0Os9jXlKjJYRCw2MzeA5YC/xUOqSEECULJQEQOVSnVYUYqVnL3ue5+LnAxwRShmwsalYhIB5SprjKOOsyoTVk/Z2Y3EzRhNeBTBY2qHaqtq2fk3OUcP3MZI+cuj6UMUURKWynVYUYtVvoYMN/dizeYeDuS3FGuvqGRy+57ibWb32XOuKExRycipaKU6jDVQ7oIRs5dnra5qwHXT/yE6itEpGjUQ7qEZKpMclBPahEpSeohXQTZKpPUk1pESlHU5KAe0odgxpjBWIafqSe1iJSi1vSQrgCeLVhE7dC4YVVccNqAlAShntQiUqqizufQCNzY/N7MXgB+Vqig2qM544ZSfexRJdEKQUSkJVGbskoelEpP6lIZ2EtESlfW5GBmX0m3GDiqMOFIoaXrczFrycsAShAisl9LTw6DMiz/ZbaNzOwO4CRgmbvPyWUdM/s58Ii7P9RCbNIKxR7YS08pIm1T1uTg7lfmukMzGw90dvcRZnanmQ1K7lmdaR0z+xTQt1QSQ3u8sBVzYC89pYi0XVFbK+WiBlgUvn4cGBVlHTMrA24D3jSz89Pt2MymmtlaM1u7ffv2vAadrJTGVc+nYg7spSlURdquQiSH7kDzFfRdoE/EdaYAvwd+Agw3s28nb+Tut7p7tbtX9+7dO++BJ2qvF7ZiDuxVSsMPi0huCpEcdgPNt6E9Mhwj3TrDCCYU2gb8J3BmAWKLrL1e2MYNq+Ka8UOpqijHgKqKcq4ZP7QgxTylNPywiOSmEE1Z1xEUJa0GTgbS3WqnW2cPwbhNANXEPGdEe54bulhNameMGXxQnQOo459IW1GI5FALPGNmlcBY4EtmNsfdL8+yzmnAh8CdZvYloAz4YgFii0wXtkNXSsMPi0huIg/ZndNOzXoB5wArw2KiVq2TTTGG7G6PrZVEpGOLOmR3QZJDMbSl+RxEREpFXudzEBGRjkVjK7WCiptEpL1TcsiRev2KSEegYqUctdfOcSIiiZQccpSu7wO0/c5xIiKJVKyUQbp6BQjGK0/XvquiWxkj5y5XPYSItAtKDmlkqlc4vKxT2sQAsPuDv/Henn0HrQ+qhxCRtknFSmlkqldovvins+/Dg9OG6iFEpC1TckgjX/UHqocQkbZKySGNTIPrVZSXpR3uuqK8LKf9iIiUOtU5JKmtq+fd9/+a9mf/ePL/ovrYo9JWVJfKIH3qoCci+aDkkOBARfSHaX/+1MbtzBmXee6DQlyUc7nYq4OeiOSLipUSXLF0Q0pFdKJMdQiFulvPdapSddATkXxRcgjV1tXT0Ji5NRKkr0Mo5FzTuV7s2+vsdSJSfEoOoZburjPVIRTybj3Xi72m5RSRfFFyCGW7u+7VrSzjPMuFvFvP9WI/Y8zgtK2pNHudiORKFdKhTHNG9+pWRt0PP5Pzdodyt95ch1Hf0JgyXEe2i72m5RSRfNFMcKHklj4AZZ2MHocfRsOefRkvtOm2Ky/rnPFJozVxNCeIKl3sReQQRZ0JTk8OoXHDqli7+V3u/e2faHLHgA+hxfGSotyt59KaKV0dRnNieG7mWfn7wCIiWSg5hGrr6rlvTZAYILggN2UYLyn5wj5uWFXe+h6oxZGIlAJVSIeufGgD+5paLmLL9SJ95UOpfSeytWZSiyMRKQUFSQ5mdoeZrTKzy3Ndx8z6mFldIeLKJtuIq4lyuUjX1tVn3G+mJNPWWxzV1tUzcu5yjp+5jJFzl+elv4eIFF/ek4OZjQc6u/sIYKCZDcpxneuAkrxNzvUina2vQ6YkM25YFdeMH0pVRTlGUNfQ2srtYitkh0ARKa5C1DnUAIvC148Do4BXo6xjZmcB7wPb0u3YzKYCUwEGDBiQz5ipKC9L20PawqZClRXlnHlib+Y9tonL7nspUjPRbEVQ2ZJMtjqMUpatQ2Bb/DwiHVkhipW6A823iu8CfaKsY2ZdgB8AMzPt2N1vdfdqd6/u3bt3HkMORlxNVtbJuH7CJ3hj7j8wY8xgFq+rz+muONvQ3+3xYqnKdJH2oxDJYTcHioV6ZDhGunVmAj9394YCxJRVbV09i9cdfJE3YOLw/gc1Vc11mIxM9QdXfO6j+Qm8xKgyXaT9KERyWEdQTARwMvBmxHXOBv7ZzFYAnzCz2wsQW1qZ+hY8tXH7/vetuStuy/UHrdHWK9NF5IBC1DnUAs+YWSUwFviSmc1x98uzrHOau9/T/EMzW+Hu3yhAbGlFufC3dpiMuOoP4pj0R8N3iLQfeU8O7r7LzGqAc4CfuPs2YH0L6+xM+nlNvuPKJsqFf8aYwWmHtTjzxPzWfeRDnJP+JCbD5gQVtQJfREpHQfo5uPt77r4oTAytXqdYohSHjBtWxRf+/uALmwP3rflTyTXVLIVJf9SsVaRtUw9potcNLPvd2ynb7mtyrnxoQ5EijaYUWg2VQoLqiK699lr69u27/1+PHj348Y9/zLRp0+jXrx9Dhw7lmWeeiTtMaQM0tlIoSt1Apt7O7+3ZF0sZfyaFGEY8V6WQoDqi733ve3zve9/b/37UqFEceeSR/PGPf+TNN9/kpZde4rzzzuOPf/wjRxxxRIyRSqnTk0OelFIRSim0GlKz1vgtX76c3r1788ADD/DNb36Tww47jOrqao455hg2bdITnGSn5BCKMiZQRXlZ2m3NSFuE8p1F62MZY6gUmtCWQoLq6G644QYuvfRSAHbs2AHAX/7yFzZv3swxxxwTZ2jSBqhYieite6743EeZcf969iUM5V3WyQ56n6h5+O9ithZqFvcQHGrWGq8tW7bw+uuvM3r0aC644AK+9a1vsXF12c9pAAAPPklEQVTjRh588EFOPPFEjjvuuLhDlBKn5EDmytMrH9oQaWKf5ik9s0kcY6iU6icKKe4E1ZHdc889TJgwAYCLL76YPn368Nxzz/Hiiy9y7733xhydtAVKDmSuJG2uaE5OEOkueMl9IDIdp1B9EDpKwpFoFi9ezG233bb//ec//3kOO+wwfvOb3zB+/PgYI5O2QnNIAyPnLs94559ues50F2I48ETRyWx/kVLyvoC0x2rtNKC1dfVcsXRDyoiyhzKPdeK+lXDannfeeYchQ4bw5z//+aDlp59+Oj/4wQ8YO3ZsTJFJKYg6h7QqpIHjjs7cgqa+oZGRc5dzee3LjJy7nONmLuOy+15KaZkE8NzMs3hj7j/w0wknZ6yMzWcTz+ankHRDjR9qn4L21IltzZo1nH/++YwZM4ZHHnnkoJ/dddddXHjhhfEEViArV65k+PDhBy17+OGH6dKlixKDRKYnB+CEWQ+nvdPPVVXC3XWmu+5MTynJTw5R7tqzPfFAMLzHG3P/oVWf5RNXPp426bT2CScuO3fu5JRTTuGWW27BzJgwYQKbNm2iV69ebN++nSFDhvCP//iPLFiwIO5QRYoi6pOD6hwgL4kBUusP0hXBpBujKbmJZ9R6iZaeNlrbp6C2rj5tYohyzFLz9ttvM2fOHM4880wA+vfvz5YtW+jVqxeXXXYZ55xzTswRipQmFSsBnc3ytq+WinOi9EGIOvREtov/ofQpaM30pqXqxBNPZMKECTQ1NbF48WL27dvHSSedxBNPPMGWLVuYOnVq3CGKlCQ9OQBfPrU//7n6rbztr6VmrS018YxaL5HuKQSgV7cyfnTeR1tdedza6U1L2Y033sjs2bO57rrr2Lt3L9OnT2fJkiXU17e9OpSo1KBADoWSAzBn3NC8Jodcn0SS/4grupWlHccp+a69UB3NMo3N1Ktb253edPr06UyYMIFTTz2Vp556ismTJzNo0KB2mxziHLZd2gclhwiqKso588TePLVxe4tPBZBah5HtDi7dH3FZJ6Oss7Gv6cB+MhUTFaKjWaZ6kR+d1/amN33ttdfYs2cPH//4x+nXrx+nnHIK69at45lnnmH+/Pns3buXxsZGdu3axZIlS+ION2+yFU0qOUgUSg6hqgx3y+laEbXU4a0q4Q6/pTu4dH/E+z50KsrL6N71sKIUCaRLXteMH9ouiiS2bt3KtGnTWLVqFY2NjaxZs4Y1a9ZQWVkJwIoVK1iwYEG7a62kUXHlUCk5hKK0IoKDi3LqGxoxgkl/Mm3T0h1cpj/WnY37eOlHn2kx7kMtV86UvK4ZP7RNNVnN5IwzzmDKlCkMGTKE7t27c8MNN+xPDO1ZKQzbLm2b+jkkaM2FtqVtjp+5jHRnuLkPQtR+D5mOnS6h5dIz+lCOL6UrH78b0j6pn0MrpCu/b774J15AO5vx5VP7M2fc0BbL/FuqXI76xJJOPsqVVfzQPmlUXDlUSg5ZZKpfaHLf37ppzrihWbff/cHfUpaXdbb9F/9D+SPOx4W9IxY/LFy4kCeffJK77rpr/7Knn36aq666iieffDLGyPJLo+LKoShIcjCzO4CTgGXuPifKOmZ2JPAroDPwPjDR3fcWIr6o0t2ZJ7r3t386KDkkFzHt2fu3tHM9dO9yWKSRXluSjwv7oTy5tDW1dfXMumEBm+67hmM+OnL/iLsvvPACX/nKVxg0aFDcIYqUjLz3kDaz8UBndx8BDDSzlL+4DOtcAMx3988A24DP5ju2KBJnhGup2Wpik9V0A9VlmnN6Z4ahKXKVj9nWSmHWuGJo/n7eWv0wFadPZM/epv0DCd58881cddVVcYcoUlIK8eRQAywKXz8OjAJebWkdd/95ws97A39O2gYzmwpMBRgwYEDeAm4WpZlqosTObi09ZSTKV5FNvsqVO0LxQ/P3c8y4Wbz/ym+AA/Uzz955J08//XTMEYqUlkIkh+5Ac7fTd4FP5rKOmY0Aern76uSN3P1W4FYIWivlMWYgtws8BMNuNItazt98Z5+voQ06woU9H5q/H0vqvb61oTFlmYgUJjnsBppvjXuQvugq7TpmdhTwM+ALBYgrq9q6+ki9n+Hg1krNMpX/J/vkgCMBWjW0gcbKab2OWPEucigKkRzWERQlrQZOBtIN8Zmyjpl1Ae4HZrn75gLElVFtXT3fuX99xp9HafOfaRC8ZM//97ts2PqXnJuglvJYOW0haXWkineRfCjEkN21wGQzmw9MADaYWXKLpeR1lgFfJyhemm1mK8xsYgFiS+vKhzbQlKZVEQSd1aJcQJIrdjNxaNVcCVGH8S62tjJjXOL3A9CtizqEiWST9ycHd99lZjXAOcBP3H0bsL6FdXYCN4f/ii5TqyIILuZRLyCJ5f8tzdKWTrYijlLtrNaWBnhr/n4W9H2LFSveK7n4REpJQfo5uPt7HGiN1Op1SkFns/3t4SF6EcqMMYO57L6X0g6d0atbGR/s+zCnIo4jy8vSPnEUosw8l2KiKEmr1IqdLrzwwpR5o2tqalixYkUs8YiUIvWQBioyXHgh6MvQXLYP0SuSxw2rYu3md1m4+q2Ugfmah76OesGsravn/b1pelp3sryXmedat9FSRW8p15WISGYaeI/gAnbpfS9lXae5rDrXQeryMWrqdxatTzvPda9uZdT9MHXk1kM5Zq4D8bU0wJsG9hMpLRp4LwfjhlW1mByyle0n/yxfxSjNF950iQGgIU1dyaHeqedat9FSR7xSrSsRkeyUHEK9Moye2qwyy5NDYrl/PotRWuqUl66+4VAriFvTHyBbR7xS7l9QanUhIqWkEE1Z26RspWvNlcVRxjLKZ5PTbHfXiSO7Rtkm6p16PsZrKuT+8qWtNMEViYueHELZBsNLbg+f7W4zn8Uo2XpdJ4/s2tI2Ue/U8z0PQKnOK9CWmuCKxEHJIZTpolpVUZ7T8Nr5LEaZMWZwxrqQTMksHz2B8z1eUymO/6S6EJHsVKwUylfxRz6LUcYNq6KivCztzzIlm44yBPehynT+SqEuRKQU6MkhlM/hr/Oxn2ZXfO6jOT8JlOKdeqnRWEsi2amfQxugVjWFofMqHVHUfg5KDiIiHUjU5KA6BxERSaHkICIiKZQcREQkhZKDiIikUHIQEZEUSg4iIpJCyUFERFIoOYiISAolBxERSdFme0ib2XZgcwF2fQywowD7LWX6zB2DPnPH0NJnPtbde7e0kzabHArFzNZG6Vrenugzdwz6zB1Dvj6zipVERCSFkoOIiKRQckh1a9wBxECfuWPQZ+4Y8vKZVecgIiIp9OQgIiIplBxERCSFkkMCM7vDzFaZ2eVxx1IMZnakmT1iZo+b2YNm1iXumIrBzPqYWV3ccRSbmf3czM6LO45CM7NeZvawma01s1vijqcYwt/pZ8LXZWb2kJk9Z2Zfa+0+lRxCZjYe6OzuI4CBZjYo7piK4AJgvrt/BtgGfDbmeIrlOqA87iCKycw+BfR194fijqUIJgMLw7b+R5hZu+7nYGa9gLuB7uGibwPr3H0k8EUzO6I1+1VyOKAGWBS+fhwYFV8oxeHuP3f3J8K3vYE/xxlPMZjZWcD7BMmwQzCzMuA24E0zOz/ueIrgHeBjZlYB9Af+FHM8hdYETAR2he9rOHAtWwm0KjkqORzQHagPX78L9IkxlqIysxFAL3dfHXcshRQWm/0AmBl3LEU2Bfg98BNguJl9O+Z4Cu1Z4FjgX4A/EPw9t1vuvsvddyYsysu1TMnhgN0cKGroQQc5N2Z2FPAzoNVlk23ITODn7t4QdyBFNgy41d23Af8JnBlzPIX2I+Bid/8xsBH4aszxFFtermUd4gIY0ToOFCWdDLwZXyjFEd5J3w/McvdCDGJYas4G/tnMVgCfMLPbY46nWF4DBoavqynMgJWlpBcw1Mw6A6cCHa0zV16uZeoEFzKznsAzwG+AscBpSY9q7Y6ZfRO4GlgfLrrZ3e+LMaSiMbMV7l4TdxzFEFZI3klQvFAGfNHd67Nv1XaZ2XDgLoKipVXA5919d7xRFV7z77SZHQs8DDwJnE5wLWvKeX9KDgeEtf7nACvDR3ARkTbHzCoJnh4ea+1NrpKDiIikUJ2DiIikUHIQEZEUSg4iIpJCyUGKyswWmNnV4esrzOyKAh5rRdL7vmZ2SB3gwviPO5R9ZNl3xvhae9zkc5Dm52WJ64Xj8qRbNsDMVpjZcjO71cws11ikbTks7gCkQ7rIzH5c7IOGLdDmFvu4UcUU3wIzewHoEQ5S9wfglDTLegPfdPc/mNkjwFDgd0WOVYpITw4Sh1cIBv0DwMy6mtm9Zva0mS1sHh02vFOdZ2aPhe/XhaPI/trMfmtmF5tZpZk9a2bPmNm/ZTuomR1nZgsS3n8lPMZKM/truKxPeIznzWxWuOz48P2TwEkt7H+hmd1lZndl2d9HzOypMO5bkrZPjC/luIlPEOGTV03Uc5DhuN8AzgX+Htju7jekW+bus939D+E2RwM7sp1rafuUHCQONwHTEt5fBLzi7qOBVzkwlMdpwCp3HxO+7wb8E/Bx4CsEvV+rCIbFGAvkNBy1u98TdoRbCfxruHgWcJ+7nw6MM7Ojge8SjEv0WaClES7PA25x9+YhG9Lt71PAy+4+ClhpZpn+DqMeN+o5SHfcW4CbgeeARjP7ToZlAJjZRGCDu2/NchxpB1SsJHHYRjDmTQ2wguCueEn4s9UEFzkIEsaShO3+x913m9lmgpEoDfgbwVg6u2n5wp3CzMYCJ7h78xweg4ERZnYhwQBmlcDxwHp3/5uZvdTCLh9PGsAw3f4eAc4ysyeA1e7+YYZ9tXTc5vFzop6DdMedAmBml7r7QU8dycvMbCBBEj07yzGkndCTg8TlemB0+HoDwVMC4f8bwtdRhjyYDlxDUBSSU4/OcJiBHwBTExZvAmaGTxRzCUa1fAv4aDhWz9AWdpscc7r9jQB+6e7nEFysT8iwr3TH3Qv0DpedEy6Leg4yHjfdUCKJy8LRA+4Fvtbeh5WRgJKDxMLd64Cnw7e3E1wEVwKDgAU57Oq/gP8HLAX2mFlVDtt+n6Ci9aGw7mEwwQX8X83sOYLinP8hKNq5HHiC4OKci3T7ex34iZmtIphDI9NAeOmO+yvgWoJin9fCZVHPQdTjpjMTGAD8LDxXo1vaQNo2DZ8hIiIp9OQgIiIplBxERCSFkoOIiKRQchARkRRKDiIikkLJQUREUvx/UYT6z5mcTrsAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from statsmodels.graphics.regressionplots import plot_leverage_resid2\n",
    "fig, ax = plt.subplots(figsize=(8,6))\n",
    "fig = plot_leverage_resid2(results, ax = ax)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Other plotting options can be found on the [Graphics page.](https://www.statsmodels.org/stable/graphics.html)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Multicollinearity\n",
    "\n",
    "Condition number:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "702.1792145490062"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.linalg.cond(results.model.exog)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Heteroskedasticity tests\n",
    "\n",
    "Breush-Pagan test:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[('Lagrange multiplier statistic', 4.893213374093957),\n",
       " ('p-value', 0.08658690502352209),\n",
       " ('f-value', 2.5037159462564373),\n",
       " ('f p-value', 0.08794028782673029)]"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "name = ['Lagrange multiplier statistic', 'p-value',\n",
    "        'f-value', 'f p-value']\n",
    "test = sms.het_breuschpagan(results.resid, results.model.exog)\n",
    "lzip(name, test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Goldfeld-Quandt test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[('F statistic', 1.1002422436378139), ('p-value', 0.38202950686925286)]"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "name = ['F statistic', 'p-value']\n",
    "test = sms.het_goldfeldquandt(results.resid, results.model.exog)\n",
    "lzip(name, test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Linearity\n",
    "\n",
    "Harvey-Collier multiplier test for Null hypothesis that the linear specification is correct:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[('t value', -1.0796490077827041), ('p value', 0.28346392475394466)]"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "name = ['t value', 'p value']\n",
    "test = sms.linear_harvey_collier(results)\n",
    "lzip(name, test)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
